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Abstract The fluid-structure interaction may occur in space 
launch vehicles, which would lead to bad performance of ve¬ 
hicles, damage equipments on vehicles, or even affect astro¬ 
nauts’ health. In this paper, analysis on dynamic behavior 
of liquid oxygen (LOX) feeding pipe system in a large scale 
launch vehicle is performed, with the effect of fluid-structure 
interaction (FSI) taken into consideration. The pipe system is 
simplified as a planar FSI model with Poisson coupling and 
junction coupling. Numerical tests on pipes between the tank 
and the pump are solved by the finite volume method. Re¬ 
sults show that restrictions weaken the interaction between 
axial and lateral vibrations. The reasonable results regarding 
frequencies and modes indicate that the FSI affects substan¬ 
tially the dynamic analysis, and thus highlight the usefulness 
of the proposed model. This study would provide a refer¬ 
ence to the pipe test, as well as facilitate further studies on 
oscillation suppression. 

Keywords Fluid-structure interaction ■ Liquid oxygen • Fi¬ 
nite volume method 

Nomenclature 

A cross-sectional area 

E Young modulus of elasticity 
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e pipe wall thickness 

/ second moment of cross-sectional area 

K liquid bulk modulus 

L length 

M bending moment 

m mass 

P pressure (cross-sectional average) 

Q lateral shear force 

R (inner pipe) radius 

t time 

u pipe displacement 

u pipe velocity 

V fluid velocity (cross-sectional average) 

x lateral coordinate (out-of-plane, vertical) 

y lateral coordinate (in-plane, horizontal) 

z axial coordinate (distance along pipe) 

A t time step (numerical grid length on f-axis) 

A z element length (numerical grid length on z-axis) 

8 rotational velocity of pipe 

v Poisson ratio 

p mass density 

cr normal stress 

Subscripts 

f fluid 
t pipe 

x lateral direction (out-of-plane, vertical) 

y lateral direction (in-plane, horizontal) 

z axial direction 
0 initial value 

Superscripts 

n number of node 
* modified value 

1 Introduction 

In flights of large scale launch vehicles with liquid rocket 
engines, axial coupled vibrations have always been ob¬ 
served [1,2], The vibrations arise when the frequency of 
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pressure oscillation in the propulsion system closely matches 
that of the vehicle structure. As the structure is difficult to 
change after design, usually our interest would focus on the 
propulsion system which includes feeding pipes and engines. 

As an important part of propulsion system study, the 
dynamic character of feeding pipe system is proposed as one 
of the key points for the analysis with the fluid-structure 
interaction (FSI) taken into consideration. Although the 
effect of FSI in liquid-filled pipes has been researched 
broadly [3,4], it is necessary for large scale launch vehicles 
design and manufacture to reinforce theoretical research on 
the FSI dynamic behavior of the feeding pipe systems. 

Skalak [5] proposed the FSI theory as an extension of 
the conventional waterhammer theory. Wiggert et al. [6] 
built an FSI model to consider Poisson coupling and junc¬ 
tion coupling, which made an advancement for the FSI the¬ 
ory. Poisson coupling relates pressures in the liquid to ax¬ 
ial stresses in the pipes through the radial vibrations of the 
pipe walls. Friction coupling is the mutual friction between 
the liquid and the axially vibrating pipe wall. For solving 
the model of FSI, the method of characteristics (MOC) was 
popular [3,6], and was improved later [7], In addition, Zhang 
et al. [8] applied the finite element method (FEM) to these 
problems, Zhang and Gay [9] presented an immersed finite 
element method (IFEM) for the solution of FSI problems, 
and Giesecke [10] solved relevant FSI mechanisms by an 
MOC-FEM hybrid method. Although successfully applied 
in the field of fluid dynamic and heat transfer, the finite vol¬ 
ume method (FVM) was not used to solve FSI problems. 
With the clear physics meaning in every term of the inte¬ 
gral equations and a well-documented numerical procedure, 
the FVM guarantees the conservation of the initial differ¬ 
ential equations, as well as the conservation of variable in 
the whole solution domain by keeping characteristic variable 
conservation in every control volume. Therefore, the FVM is 
adopted in present paper to solve the FSI problems for LOX 
feeding pipe systems. 

The control equations of the FSI model were solved by 
the FVM in the time domain under the well-posed condition. 
Through numerical calculation, the dynamic response of the 
LOX pipes was studied, and the FSI mechanisms was dis¬ 
cussed. 

2 Description of mathematical model 

As shown in Fig. 1, one liquid-filled pipe element is assumed 
to be a planar model, which can vibrate in the axial (longi¬ 
tude) and lateral directions, neglecting vibration in the verti¬ 
cal direction. In other words, the fluid-filled pipe system can 
only move in the y-z, plane without consideration of motion 
in ^-direction. The parameters interested include pipe wall 
axial velocity ii z , pipe wall lateral velocity u y , fluid pressure 
P, fluid velocity V, pipe wall axial stress cr,, pipe wall lateral 
shear force Q y , pipe angular velocity 0 X and bending moment 
M x . 



Fig. 1 Element of liquid-filled pipe hint 


For simplification, some assumptions are made, in¬ 
cluding: (1) The fluid is considered as ID axial flow, with 
changes in the radial direction neglected; (2) The hoop stress 
of the pipe wall section is equally distributed; (3) The pipe 
wall is a linear elastic body. 

According to the FSI mechanisms including Poisson 
coupling and friction coupling, with the Timoshenko beam 
mechanisms taken into consideration, the motion differential 
equations were built in the axial and lateral directions for a 
liquid-filled pipe, as set up by Wiggert [6]. 

The four axial motion differential equations are written 
as Eqs. (l)-(4). 

The fluid axial motion equation with the viscosity taken 
into consideration is 


dV 

dt 


1 dP 

+ — — +R f V = 0, 
Pf oz 


(1) 


where pf is fluid density, Rf is fluid viscosity damp coeffi¬ 
cient, t is time, z is distance of ^-direction. 

The fluid continuity equation with the compressibility 
taken into consideration is 


dV 1 dP ^ <3w z 
~dz + K*~dt~ " V ~dz 


( 2 ) 


where v is Poisson ratio, K* is fluid bulk modulus correction. 
The pipe axial motion equation is 


dip 

dt 


1 dcr- 
Pt oz 


( 3 ) 


where p t is density of pipe wall, R t is pipe damp coefficient. 

The pipe stress-strain relation with the fluid pressure 
effect taken into consideration is 


Hit- 1 dcr. vR dP 

dz E dt Ee dt ’ 


where E is Young modulus, R is inner radius, e is wall thick¬ 
ness. 

The lateral motion differential equations are also writ¬ 
ten as 


dUy 1 dQy 

dt pA dz 

diiy 1 dQ y 

~dz + GA t ~dT 


( 5 ) 

( 6 ) 
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where inertia per length pA = />, A , + pfAf, and G is the pipe 


wall shear modulus. 


d0 x 1 8M X _ 1 

dt pi dz pi " 

(7) 

+ 

C»3 ^ 

Cu 

II 

O 

(8) 


where I is second moment of area, equivalent moment of the 
inertia per length pi = p t / t + p f / f . 


3 Numerical solutions by FVM 


According to the liquid-filled pipe axial motion differential 
equations (l)-(4), the axial motion difference equations can 
be achieved by time discrete integral method in a staggered 
grid. 

Take Eq. (1) as an example, the integral form in the 
control volume is written as 




dvdr = 0. 


(9) 


The selected control volume is named as control vol¬ 
ume V, as shown in Fig. 2. 



Fig. 2 Control volume V and nodes distribution 


The left first integral term in Eq. (9) can be written as 

J f* /■» c \tt r* y w+ 1 y/2 

-r-dtdv = - -dtdv = (V" +1 - V")AV. 

AvJAr dt J AV At 

For ID grid, AV - A u Az, the above equation may be 
written as 

y/l+t _ yn 


r r dV C 

—dvdr= 

JAV J At dt J A 


AV 
- /t/M+l 


At 


-dfdv 


= (V’’ +l - V’})A U Az. 


( 10 ) 


In Eq. (10) Aj j is the area of the control volume inter¬ 
face of the grid. For ID, it is equal to the section area of the 
fluid flowed in the pipe A f . Az is the ID grid width 6z, called 
as the space step. 

In the second integral term in Eq. (9), the value of the 
left and right interfaces of the control volume are the nodes 


£) Springer 


value of the P control volume, written as Pj and P/- \. 

P - P 

1 p *- u; 


r r l dp ir 
dTdt= - 

JAt JAV Pf OZ Pf JAi 


-A/.ydzdr 


At Az 

- r p ' ~ Pi 1 Aijdzdt. 
Pf J At Az 


Pf 

Using the method of weighed average and integration 
we can obtain a relation as 


n _L _8P_ 

v Pf dz 


dVdt = ^^[0(P” +1 - P n ! + _\) 
Pf 


+(1 - ff)(P'l - PIJ]. 


(ID 


In the third term in Eq. (9), V, at the node is used in¬ 
stead of the average of the element the node belongs to. 

n RfVdvdt = R f A u AzAt\eV'; +l + (1 - 0)V?]. (12) 

v 

Substituting Eqs. (10)—(12) into Eq. (9), the differenc¬ 
ing scheme of Eq. (1) can be written as 

.At , , 

(p f + p f R f At0)V n+l + —6»(F, !+1 - F/+, 1 ) = MV". P"). (13a) 
Az 

The right side term fi(V n , P") is a function of the value 
of one node at the n time layer, which can be written as 


Mv n ,p n ) = [p f +p i R i At{e-\)]v? 

+^(1 -OXPI, - P']). (13b) 

Equation (2) is treated with the same discrete method, 
but in control volume P as shown in Fig. 3. 



Fig. 3 Control volume P and nodes distribution 


According to the basic principle for discretizing the 
above equations, the axial motion difference equations can 
be obtained as 
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(Pt+PfRfmv'^ 1 + —8(P‘; +l - p n +{) 


= f\(V n ,P n ), 
w, 1 - V " +1 ) + 


pn+l 

K* At 1 


n+l 
i +1 


2 ve[(ii z y. 

-(uz)T 1 ]=f2(v n ,P n ,u n z ), 

- Pt R t AzAt8V? +l +p l Az(u z )- + ' - Ate\(cr z ) 

—A u AzP'j +1 + A u EAte\(u z )^l - 0 u z )1 +l ] 


n+l 

I 


(14a) 


-A u Az(o- z )l +1 =f 4 {P n ,u n z ,(f' z ). 

The related functions at one node at the n time layer can 
be written as 

Mv n ,p n ) = \p f +p f RfAt(e- i)]yf 
MV", P\ u’l) = {6- 1)(V? +1 - V?) + -^P’l 

+ 2y(\-0)\(u z )i,-(u 7 y;\- (14b) 

MV", it", M) - p t R t AzAt(l - 0)vp + p t Az(uj; 

+Af(i - 0)f(cr z y; - (cr z )?_, \, 

f 4 (P", it", (X’l) = —A/ jAzP" + A, jEAt(9 - 1) 
e 

x[(u z y; +l -(u z y;\-A u A Z (cr z r r 


i 


GA t At 


(Qy) n , + 


(0-1X0*)?, 


0?, M'l) = ^%)? + ^j(0 - 1) (15b) 

x[(M x y; - (m x )?_j] + ^(1 - 0)(Qy)", 

= {6- 1)[ (0 X )" +1 - ( 0 x y;\ + 

To solve the partial differential equations needs some 
conditions, called as solvability conditions. Usually, in re¬ 
gion Q (x e Q), to satisfy the equation, the relationship of u 
at the boundary of region Q is boundary condition. 

The control equation system consists of the above eight 
partial differential equations. The conditions for constitut¬ 
ing a well-posed problem are the initial conditions and the 
boundary conditions in different forms, which satisfy con¬ 
sistent, astringency and stability requirements. The initial 
conditions do not affect the solution, whereas most bound¬ 
ary conditions affect the form of discrete equations and the 
calculation method. 

In the staggered grid, one-extra-layer nodes are set at 
the outside of the boundary, as shown in Fig. 4. Calculation 
is developed from the inner node {1-2), the boundary node 
is used only for giving or calculating the boundary value. 

(1) Physical boundary coincides with the control vol¬ 
ume boundary. 

(2) Nodes of outside layers are used to define, calculate 
or store the boundary values. 


In the same manner, the lateral motion difference equa¬ 
tions can be written as 

- (QyT-i] = q> 

e 


A z 


pA 

[«/ - w 1 ] 


i 


GA t At 


{Qy)T l + 0(0*)r 


= f 6 (u;,Q';,o ! D, 


A ute, b 


At 


( 0 ,r' + - (M,y;n] 


(15a) 


AfAi 

Pi 


8(Qy)i = MQ", 0?, M'l), 


The related functions at one node at the n time layer can 
be written as 

m,Q;) = ^(u y y; + ^(e-\) 

x[(Qy)1 ~ (Qy)Ul 

/ 6 («?,e?,^) = ^[(« y )? +I -K); ! ] 




/= I /= 1 1=2 i=2 ■■■ i=n -1 I=n i=n 


Fig. 4 Boundary gird configuration 


These works minimize the changes of the discrete equa¬ 
tions and the integral for the boundary of the control volume 
by the boundary condition. 

The number of the physical and numerical boundaries 
at every boundary node is equal to the number of variables, 
to ensure the solution of these nodes. With consideriation 
of the flow direction, the number of the physical boundary 
may be equal to the number of boundary node with input 
character. These two boundaries are shown in algebra form. 
If the number of the physical boundary is greater than that 
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for a well-posed problem, the extra ones are considered as 
additive boundary conditions. 

The calculation form of the numerical boundary condi¬ 
tion may be the same as the discrete form of the inner nodes. 
Due to the staggering of the discrete form of variables in 
the staggered grid system, the physical boundaries can be 
matched to the nodes of the control volume V, thus, the cal¬ 
culation form of the numerical boundaries is the same as the 
form of nodes of control volume P. The boundaries of the 
discrete form are self-contained when the number of physi¬ 
cal quantities and the number of variables are equal. 


4 Numerical example for testing 


The experiment of University of Dundee is simulated by the 
method proposed in the present paper. The experimental 
apparatus consists of one long and one short stainless steel 
pipes screwed onto a rigid 90° elbow [3]. Inner diameter of 
the pipes is 52 mm, and the wall thickness is 3.9 mm. The 
pipes is filled with pressurized tap water (initial static pres¬ 
sure 2 MPa), with both ends closed. The pipes can move 
freely in a nearly horizontal plane because it is suspended on 
three long steel wires. Vibrations are generated by striking 
the closed end of the long pipe axially with a solid steel rod 
of 5 m length. Just before the impact, the rod moves with a 
constant velocity of 0.809 m/s. 

The calculated and measured pressure results close to 
the elbow are shown in Fig. 5. The agreement between the 
numerical simulation of the FSI model by FVM and the ex¬ 
perimental result is good. The results verified the rationality 
of the simulation of the FSI model by FVM, and the validity 
of the simulation code. These would facilitate further simu- 
lational studies. 



Fig. 5 Simulated and measured pressures close to the elbow 


5 Numerical simulation of LOX feeding pipe before 
pump 


For a typical launch vehicle, the FSI usually occurs in the 
£) Springer 


LOX feeding pipe systems. The system of LOX feeding 
pipe before pump is now analyzed. As shown in Fig. 6, the 
right Pipe 3 is the model of LOX feeding pipe. The middle 
Pipe 2 is the model whose two middle restrictions are absent. 
The left Pipe 1 is a straight pipe just for comparison, whose 
length is the same as the length of the other two pipes. 
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Fig. 6 Feeding Pipes 1, 2 and 3 


The inlet of the pipes is also the outlet of the storage. 
The storage pressure at the outlet is Pq = 0.4 MPa and the 
volume flux is Qq = 0.2 m 3 /s. The boundary conditions are 


m,(0, t) = 0, P(0,t) = P 0 , 

0,(0, t) = 0, M x (0, t) - 0. 


The end of the pipes is fixed. A spring pulse is set at 
the end, which is 1% flux reduction during 0.01 s. The flux 
oscillation is expressed as 


0(0 = 00 —y^U-T), 



t c = 0.01 s, 


E m = 1.5. 


The initial stress, shearing force and bending moments 
are analyzed by the statics balance condition. 

The pressure response to spring disturbance at the end 
of pipe is shown in Fig. 7. 

For the straight Pipe 1, since the FSI effect is consid¬ 
ered, in pressure response Curve 1 there appear many hack¬ 
les oscillation compared with the traditional water hammer 
model. These hackles reflect the effect of the stress wave. 

For Pipe 2, because of the junction coupling, the axial 
vibrations interact with the lateral vibrations, and thus Curve 
2 becomes very complex. The spring pulse is applied in the 
axial direction, however, due to the axial and lateral balance, 
the vibration form is changed to the axial and lateral mixed 
vibration. 
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Fig. 7 Curves of the pressure response to spring disturbance at the 
end of pipe 

For adding hinge bearings at the elbow as in Pipe 3, the 
effect of the interaction between the axial and lateral vibra¬ 
tions is reduced, which means less elbow junction coupling. 
The curve of Pipe 3 is similar to that of Pipe 1. 

All power spectrum density estimates of these curves 
are shown in Figs. 8-10. The base frequency of Curve 1 is 
9.26 Hz, and that of Curve 3 is 9.26 Hz, too. However, the 
base frequency of Curve 2 is 8.23 Hz, as shown in Fig. 9. Be¬ 
cause the results of Curve 2 are affected by Poisson coupling 
and junction coupling, the pressure oscillation peak comes 
later than in other two curves. In addition, there are some 
high frequencies components in the power spectrum density 
estimate shown in Fig. 9, which reflect the effects of stress 
wave, shear wave and bending wave. 

In addition, the power at the base frequency is also 
reduced for Curve 2. As the axial and lateral vibrations are 
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Fig. 8 Power spectrum density estimate of Curve 1 
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Fig. 9 Power spectrum density estimate of Curve 2 
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Fig. 10 Power spectrum density estimate of Curve 3 

coupled, the amplitude power at the frequency of lateral vi¬ 
bration is increased, and thus that of axial low frequency vi¬ 
bration is reduced. 

6 Summary 

The simulation results of the FSI model by FVM are found 
in close agreement with the experiment results reported in 
literature, which validats the code of FVM. In the present 
simulation of LOX feeding pipe system, the axial and lat¬ 
eral vibrations are found to have severe interaction for the 
case without any restraint at the elbow. The simulation re¬ 
sults show a complex pressure response curve resulting from 
the effect of elbow conjunction coupling. However, if the 
elbow vibration is restricted, the pressure response curve is 
similar to that of a straight pipe, due to the suppression ef¬ 
fect of the elbow conjunction coupling. For complex feeding 
pipe systems used in space launch vehicles, if the dynamic 
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characteristics of pipes are to be analyzed accurately, the dis¬ 
tribution of pipes, the elbow pipes and the restrictions are all 
needed to consider in detail. 
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